Taller INLA covid
1 Analisis Predictivo
El análisis predictivo emplea datos históricos para predecir eventos futuros. Normalmente, los datos históricos se utilizan para crear un modelo que capture diversos patrones temporales. Este modelo predictivo se usa entonces con los datos actuales para predecir lo que acontecera. Existen diversas metodologias para llevar acabo predicciones, en este taller optamos por el enfoque bayesiano para llevar acabo la prediccion.
1.1 Enfoque bayesiano : Distribucion predictiva aposteriori
El enfoque bayesiano, en el contexto del analisis predictivo, enfoca su analisis en el calculo de la distribucion predictiva aposteriori. Esta distribucion emplea la informacion de los datos asi como la distribucion aposteriori de los parametros de interes . Empleando estos insumos la distribucion predictiva podemos generar datos futuros; es decir, datos faltantes en un determinando horizonte.
\[\begin{align*} f(y_{pred}|y_{obs}) = \int_{\theta} f(y_{pred},\theta|y_{obs})d\theta =\int_{\theta} f(y_{pred}| \theta y_{obs})f(\theta|y_{obs})d\theta \end{align*}\]
Actualmente existen diversas metodologias para llevar acabo el calculo numerico de esta distribucion las cuales en R corresponden al empleo de un determinado paquete en. Algunos paquetes ampliamente conocidos son los siguientes :
- STAN
- Winbugs
- Jags
Relativamente hace no mucho tiempo la metolodogia de la Aproximación anidada integrada de Laplace o INLA(por sus siglas en ingles) ha cobrado importancia por su relativa eficiencia para estimar modelos complejos.
1.2 Prediccion como datos faltantes
El analisis de prediccion como en la practica puede ser entendida como predecir datos desconocidos en un determinado horizonte. En ese sentido estos datos a futuro pueden ser interpretados como datos faltantes o missings values, los cuales son el objeto de la predicción.
| Año | Semana | Numero de muertes |
|---|---|---|
| 2019 | 1 | 0 |
| 2019 | 2 | 1 |
| 2019 | 3 | 4 |
| 2020 | 1 | 0 |
| 2020 | 2 | 1 |
| 2020 | 3 | 4 |
| 2021 | 1 | NA |
| 2021 | 2 | NA |
| 2021 | 3 | NA |
Para llevar acabo predicciones empleando INLA, una aproximacion es justamente la anterior. Una vez dispuestos los datos como en la tabla anterior, basta con ajustar un modelo tipico en INLA, ya que INLA internamente calculara la distribucion predictiva aposteriori .
2 Analisis Descriptivo
library(tidyverse)
library(INLA)
library(yardstick)
library(gt)
library(innovar)
library(spdep)
library(reshape2)
db <- readRDS("db_excess_proc.rds")
db.lima <- db %>%
filter(prov=="LIMA")
db.lima %>%
ggplot()+
geom_line(aes(x=date,y=n),color="#011f4b")+
geom_point(aes(x=date,y=n),color="black")+
facet_wrap(vars(sexo))+
ylab("Numero de muertes")+
xlab("")+
ggtitle("Evolucion del Numero de muertes por género")+
theme_bw()+
theme(plot.title = element_text(hjust = 0.5,size = "22"),
strip.background =element_rect(fill="#011f4b"),
strip.text = element_text(color = "white",face = "bold",size = "5pts")) 
2.1 Modelamiento con INLA
Para llevar acabo el ajuste de modelos empleamos la funcion inla(), la cual tiene diversos parametros, algunos de los mas importantes son :
data: Un objeto tipicamente de la clasedata.frameel cual representa los datos para ajustar el modeloformula: Un objeto de la clase formula que especifica la ecuacion que pretendemos ajustar como por ejemploinla(target ~ 1 + x1). En la formula podemos especificar efectos lineales (introduciendo el nombre de lavariable) o no lineales (empleandof()).verbose: Una variable del tipoboolean, la cual señala si se desea mostrar en consola el proceso de convergencia
db.femenino <- db %>%
filter(sexo=="FEMENINO") %>%
mutate(n=ifelse(annus == 2019,NA,n))
db.masculino <- db %>%
filter(sexo=="MASCULINO") %>%
mutate(n=ifelse(annus == 2019,NA,n))
db.n.femenino <- db.lima %>%
filter(sexo=="FEMENINO") %>%
mutate(n=ifelse(annus == 2019,NA,n))
db.n.masculino <- db.lima %>%
filter(sexo=="MASCULINO") %>%
mutate(n=ifelse(annus == 2019,NA,n))lima.m.m0 <- inla(n ~ 1 + annus,
verbose = F,
data = db.n.masculino
)
lima.f.m0 <- inla(n ~ 1 + annus,
verbose = F,
data = db.n.femenino
)Los parámetros antes detallados son los esenciales para ejecutar el ajuste de un modelo empleando INLA.Sin embargo, algunos parametros extra a tener en consideracion son los siguientes :
family:Objeto de clasecharacter.Este parametro es crucial, pues determina la distribucion de la variable objetivo, por defecto se encuentra enfamily=Gaussian.control.compute:Objeto de la claselistpermite especificar el calculo de criterios de informacion tales comoaic,dic,waic.
lima.m.m0 <- inla(n ~ 1 + annus,
verbose = F,
data = db.n.masculino,
control.compute = list(dic=TRUE, waic=TRUE,cpo=TRUE),
control.predictor = list(link = 1),
)
lima.f.m0 <- inla(n ~ 1 + annus,
verbose = F,
data = db.n.femenino,
control.compute = list(dic=TRUE, waic=TRUE,cpo=TRUE),
control.predictor = list(link = 1),
)2.2 Prediccion: Enfoque Temporal
En el analisis de prediccion (forecast) tipicamente se emplea datos en el formato de series de tiempo para llevarla acabo. Para ello comunmente se asume que las series de tiempo pueden ser descompuestas del siguiente modo :
\[\begin{align*} y_{t} = S_{t}+T_{t}+e_{t} \end{align*}\]
donde :
\(S_{t}\) : Estacionalidad
\(T_{t}\) : Tendencia
\(e_{t}\) : Ruido blanco
2.2.1 Ajuste y prediccion
Los componentes antes detallados pueden ser modelados en el contexto de INLA empleando distribuciones apriori para cada componente tales como los siguientes:
- linear :
variable
- AR1 (Proceso autoregresivo de 1er orden) :
f(variable, model = "ar1")
- RW1 (Paseo aleatorio de 1er orden) :
f(variable, model = "rw1")
- RW2 (Paseo aleatorio de 2do orden) :
f(variable, model = "rw2")
Dado el comportamiento que exibe la series que prentedemos modelar, para los componentes :
- Tendencia : Emplearemos una distribucion apriori del tipo autoregresivo(AR1) o linear
- Estacionalidad : Emplearemos una distribucion apriori de tipo de paseo aleatorio de 1er o 2do orden
2.2.1.1 Hombres
#inla.list.models("likelihood")
lima.m.m1 <- inla(n ~ 1 + f(annus,model = "ar1"),
control.compute = list(dic=TRUE, waic=TRUE,cpo=TRUE),
control.predictor = list(link = 1),
verbose = F,
family = "nbinomial",
data = db.n.masculino
)
lima.m.m2 <- inla(n ~ 1 + f(week,model = "rw1"),
control.compute = list(dic=TRUE, waic=TRUE,cpo=TRUE),
control.predictor = list( link = 1),
verbose = F,
family = "nbinomial",
data = db.n.masculino
)
lima.m.m3 <- inla(n ~ 1 + f(annus,model = "ar1") + f(week,model = "rw1"),
control.compute = list(dic=TRUE, waic=TRUE,cpo=TRUE),
control.predictor = list( link = 1),
verbose = F,
family = "nbinomial",
data = db.n.masculino
)fit.m.m1 <- lima.m.m1$summary.fitted.values$mean
fit.m.m2 <- lima.m.m2$summary.fitted.values$mean
fit.m.m3 <- lima.m.m3$summary.fitted.values$mean
n.m <- db.n.masculino$n
datos <- list("AR(1)"=fit.m.m1,"RW(1)"=fit.m.m2,"AR(1)RW(1)"=fit.m.m3) %>%
as.data.frame() %>%
melt( variable.name = "modelo",
value.name = "fit") %>%
group_by(modelo) %>%
mutate(actual = n.m,
date = db.n.masculino$date) %>%
ungroup() %>%
mutate(min_inter = c(lima.m.m1$summary.fitted.values$`0.025quant`,
lima.m.m2$summary.fitted.values$`0.025quant`,
lima.m.m3$summary.fitted.values$`0.025quant`),
max_inter = c(lima.m.m1$summary.fitted.values$`0.975quant`,
lima.m.m2$summary.fitted.values$`0.975quant`,
lima.m.m3$summary.fitted.values$`0.975quant`)
)
datos %>%
ggplot() +
geom_line(aes(x=date,y=actual)) +
geom_line(aes(x=date,y=fit),color="#011f4b") +
geom_ribbon(aes(x=date,y=fit,ymin=min_inter, ymax=max_inter),
alpha=0.4, #transparency
fill="#011f4b") +
facet_wrap(vars(modelo)) +
theme_bw()+
xlab(" ") +
ylab("Numero de casos")+
geom_vline(xintercept = as.Date("2019-01-07"), linetype="dotted",
color = "black", size=0.7)+
theme(plot.title = element_text(hjust = 0.5,size = "20pts"),
strip.background =element_rect(fill="#011f4b"),
strip.text = element_text(color = "white",face = "bold",size = "5pts")) 
2.2.1.2 Mujeres
#inla.list.models("likelihood")
lima.f.m1 <- inla(n ~ 1 + f(annus,model = "ar1"),
control.compute = list(dic=TRUE, waic=TRUE,cpo=TRUE),
control.predictor = list( link = 1),
verbose = F,
family = "nbinomial",
data = db.n.femenino
)
lima.f.m2 <- inla(n ~ 1 + f(week,model = "rw1"),
control.compute = list(dic=TRUE, waic=TRUE,cpo=TRUE),
control.predictor = list(link = 1),
verbose = F,
family = "nbinomial",
data = db.n.femenino
)
lima.f.m3 <- inla(n ~ 1 + f(annus,model = "ar1") + f(week,model = "rw1"),
control.compute = list(dic=TRUE, waic=TRUE,cpo=TRUE),
control.predictor = list( link = 1),
verbose = F,
family = "nbinomial",
data = db.n.femenino
)fit.f.m1 <- lima.f.m1$summary.fitted.values$mean
fit.f.m2 <- lima.f.m2$summary.fitted.values$mean
fit.f.m3 <- lima.f.m3$summary.fitted.values$mean
n.f <- db.n.femenino$n
datos <- list("AR(1)"=fit.f.m1,"RW(1)"=fit.f.m2,"AR(1)RW(1)"=fit.f.m3) %>%
as.data.frame() %>%
melt( variable.name = "modelo",
value.name = "fit") %>%
group_by(modelo) %>%
mutate(actual = n.m,
date = db.n.femenino$date) %>%
ungroup() %>%
mutate(min_inter = c(lima.f.m1$summary.fitted.values$`0.025quant`,
lima.f.m2$summary.fitted.values$`0.025quant`,
lima.f.m3$summary.fitted.values$`0.025quant`),
max_inter = c(lima.f.m1$summary.fitted.values$`0.975quant`,
lima.f.m2$summary.fitted.values$`0.975quant`,
lima.f.m3$summary.fitted.values$`0.975quant`)
)
datos %>%
ggplot() +
geom_line(aes(x=date,y=actual)) +
geom_line(aes(x=date,y=fit),color="#011f4b") +
geom_ribbon(aes(x=date,y=fit,ymin=min_inter, ymax=max_inter),
alpha=0.4, #transparency
fill="#011f4b") +
facet_wrap(vars(modelo)) +
theme_bw()+
xlab(" ") +
ylab("Numero de casos")+
geom_vline(xintercept = as.Date("2019-01-07"), linetype="dotted",
color = "black", size=0.7)+
theme(plot.title = element_text(hjust = 0.5,size = "20pts"),
strip.background =element_rect(fill="#011f4b"),
strip.text = element_text(color = "white",face = "bold",size = "5pts")) 
2.2.2 Performance
2.2.2.1 Hombres
perform.metrics <- metric_set(mae, rmse,msd)
tbl.yrd.m <- datos %>%
group_by(modelo) %>%
perform.metrics(truth = n.m,
estimate = fit)
tbl.yrd.m %>%
pivot_wider(id_cols = modelo,
names_from = .metric,
values_from = .estimate) %>%
gt() %>%
tab_header(
title = md("Métricas de ajuste"),
subtitle = "Hombres") %>%
data_color(
columns = vars(rmse,mae,msd),
colors = scales::col_numeric(
palette = c(
"#011f4b","white"),
domain = NULL)
)| Métricas de ajuste | |||
|---|---|---|---|
| Hombres | |||
| modelo | mae | rmse | msd |
| AR.1. | 40.66904 | 50.09989 | 27.70859 |
| RW.1. | 34.64344 | 42.68089 | 27.66907 |
| AR.1.RW.1. | 30.38785 | 37.24076 | 27.67888 |
2.2.2.2 Mujeres
perform.metrics <- metric_set(mae, rmse,msd)
tbl.yrd.f <- datos %>%
group_by(modelo) %>%
perform.metrics(truth = n.f,
estimate = fit)
tbl.yrd.f %>%
pivot_wider(id_cols = modelo,
names_from = .metric,
values_from = .estimate) %>%
gt() %>%
tab_header(
title = md("Métricas de ajuste"),
subtitle = "Mujeres") %>%
data_color(
columns = vars(rmse,mae,msd),
colors = scales::col_numeric(
palette = c(
"#011f4b","white"),
domain = NULL)
)| Métricas de ajuste | |||
|---|---|---|---|
| Mujeres | |||
| modelo | mae | rmse | msd |
| AR.1. | 30.25786 | 35.90613 | 0.10282570 |
| RW.1. | 24.32470 | 30.10790 | 0.06330288 |
| AR.1.RW.1. | 18.04854 | 22.13154 | 0.07311218 |
2.3 Prediccion: Enfoque Espacio-Temporal
El enfoque espacio-temporal, implica complementar el enfoque anterior complejizando el modelo tomando en cuenta la dimesion espacial. Para ello requerimos incluir efectos aleatorios espaciales para controlar por tal dimension de analisis. Cabe mencionar que modelos espaciales areales empleando INLA, requieren necesariamente la creacion de un grafo(graph) o una estructura de datos espaciales que captura la disposicion de las areas como una matriz de pesos espaciales.
\[\begin{align*} y_{t} = S_{t}+T_{t}+\nu_{t}+u_{t}+e_{t} \end{align*}\]
donde:
\(\nu_{t}\) : efectos estructurados
\(u_{t}\) : efectos no estructurados
Los 2 nuevos componentes \(\nu_{t},u_{t}\) pretenden tomar en consideracion en la prediccion la correlacion espacial presente en los datos . Al respecto :
Los efectos no estructurados corresponden a efectos aleatorios que pretenden controlar caracteristicas no observadas de cada area bajo estudio. Para modelarlas podemos emplear
f(area, model = "iid")Los efectos estructurados son efectos aleatorios que toman en cuenta explicitamente la estructura espacial
- Pueden ser modelados en INLA de diversas maneras:
- Empleando efectos espaciales besag :
f(area, model = "besag")
- Empleando efectos espaciales propios de besag :
f(area, model = "besagproper")
- Empleando efectos espaciales besag :
- Pueden ser modelados en INLA de diversas maneras:
2.3.1 Preprocesamiento
2.3.2 Ajuste y prediccion
2.3.2.1 Hombres
peru.m.m5 <- inla(n ~ 1 + f(annus,model = "linear") +
f(week,model = "rw1") +
f(id.sp, model = "besag", graph = W.peru),
data = db.m.sp,
control.compute = list(dic = TRUE,
waic = TRUE,
cpo = TRUE),
control.predictor = list( link = 1)
)
2.3.2.2 Mujeres
peru.f.m5 <- inla(n ~ 1 + f(annus,model = "ar1") +
f(week,model = "rw1") +
f(id.sp, model = "besag", graph = W.peru),
data = db.f.sp,
control.compute = list(dic = TRUE,
waic = TRUE,
cpo = TRUE),
control.predictor = list( link = 1)
)
2.3.3 Performance
2.3.3.1 Hombres
perform.metrics <- metric_set(mae, rmse,msd)
tbl.yrd.m <- datos %>%
group_by(modelo) %>%
perform.metrics(truth = n.m,
estimate = fit)
tbl.yrd.m %>%
pivot_wider(id_cols = modelo,
names_from = .metric,
values_from = .estimate) %>%
gt() %>%
tab_header(
title = md("Métricas de ajuste"),
subtitle = "Hombres") %>%
data_color(
columns = vars(rmse,mae,msd),
colors = scales::col_numeric(
palette = c(
"#011f4b","white"),
domain = NULL)
)| Métricas de ajuste | |||
|---|---|---|---|
| Hombres | |||
| modelo | mae | rmse | msd |
| AR.1. | 30.25786 | 35.90613 | 0.10282570 |
| RW.1. | 24.32470 | 30.10790 | 0.06330288 |
| AR.1.RW.1. | 18.04854 | 22.13154 | 0.07311218 |
2.3.3.2 Mujeres
perform.metrics <- metric_set(mae, rmse,msd)
tbl.yrd.m <- datos %>%
group_by(modelo) %>%
perform.metrics(truth = n.m,
estimate = fit)
tbl.yrd.m %>%
pivot_wider(id_cols = modelo,
names_from = .metric,
values_from = .estimate) %>%
gt() %>%
tab_header(
title = md("Métricas de ajuste"),
subtitle = "Mujeres") %>%
data_color(
columns = vars(rmse,mae,msd),
colors = scales::col_numeric(
palette = c(
"#011f4b","white"),
domain = NULL)
)| Métricas de ajuste | |||
|---|---|---|---|
| Mujeres | |||
| modelo | mae | rmse | msd |
| AR.1. | 30.25786 | 35.90613 | 0.10282570 |
| RW.1. | 24.32470 | 30.10790 | 0.06330288 |
| AR.1.RW.1. | 18.04854 | 22.13154 | 0.07311218 |